library(data.table)
library(magrittr)
library(here)
dir = here()
# model
wd = gsub("Analysis", "Hierarchy Model/results", dir)
setwd(wd)
mod = fread("inequality_energy.csv")
# fit beta
fit_beta = function(energy, gini){
error = abs(log(energy) - log(mod$energy)) + abs(gini - mod$pay_gini)
fit.id = which.min(error)
beta = mod$pay_exponent[fit.id]
return(beta)
}
# Modern
################################################################3
wd = gsub("Analysis", "Empirical Data/Energy Inequality Merge", dir)
setwd(wd)
modern = fread("energy_gini.csv")
modern$beta = mapply(fit_beta, energy_gini$energy, energy_gini$`gini_50%`)
modern_top = fread("energy_top1.csv")
wd = gsub("Analysis", "Empirical Data/Energy Inequality Merge", dir)
setwd(wd)
modern = fread("energy_gini.csv")
modern$beta = mapply(fit_beta, energy_gini$energy, energy_gini$`gini_50%`)
modern$beta = mapply(fit_beta, modern$energy, energy_gini$`gini_50%`)
modern = fread("energy_gini.csv")
modern$beta = mapply(fit_beta, modern$energy, modern$`gini_50%`)
modern_top = fread("energy_top1.csv")
fit_beta_top = function(energy, top){
error = abs(log(energy) - log(mod$energy)) + abs(gini - mod$pay_top_1)
fit.id = which.min(error)
beta = mod$pay_exponent[fit.id]
return(beta)
}
modern_top = fread("energy_top1.csv")
modern_top$beta = mapply(fit_beta_top, modern_top$energy_pc, modern_top)
error = abs(log(energy) - log(mod$energy)) + abs(top - mod$pay_top_1)
fit.id = which.min(error)
fit_beta_top = function(energy, top){
error = abs(log(energy) - log(mod$energy)) + abs(top - mod$pay_top_1)
fit.id = which.min(error)
beta = mod$pay_exponent[fit.id]
return(beta)
}
modern_top$beta = mapply(fit_beta_top, modern_top$energy_pc, modern_top)
fit_beta_top = function(energy, top){
error = abs(log(energy) - log(mod$energy)) + abs(top - mod$pay_top_1)
fit.id = which.min(error)
beta = mod$pay_exponent[fit.id]
return(beta)
}
wd = gsub("Analysis", "Empirical Data/Energy Inequality Merge", dir)
setwd(wd)
modern = fread("energy_gini.csv")
modern_top = fread("energy_top1.csv")
modern_top$beta = mapply(fit_beta_top, modern_top$energy_pc, modern_top)
modern_top$beta = mapply(fit_beta_top, modern_top$energy_pc, modern_top$`top_1_50%`)
plot(modern_top$energy_pc, modern_top$beta)
plot(modern_top$energy_pc, modern_top$beta, log = "x")
library(data.table)
library(magrittr)
library(here)
dir = here()
# model
wd = gsub("Analysis", "Hierarchy Model/results", dir)
setwd(wd)
mod = fread("inequality_energy.csv")
# fit beta
fit_beta = function(energy, gini){
error = abs(log(energy) - log(mod$energy)) + abs(gini - mod$pay_gini)
fit.id = which.min(error)
beta = mod$pay_exponent[fit.id]
return(beta)
}
# fit beta top
fit_beta_top = function(energy, top){
error = abs(log(energy) - log(mod$energy)) + abs(top - mod$pay_top_1)
fit.id = which.min(error)
beta = mod$pay_exponent[fit.id]
return(beta)
}
# Modern
################################################################3
wd = gsub("Analysis", "Empirical Data/Energy Inequality Merge", dir)
setwd(wd)
modern = fread("energy_gini.csv")
modern$beta = mapply(fit_beta, modern$energy, modern$`gini_50%`)
modern_top = fread("energy_top1.csv")
modern_top$beta = mapply(fit_beta_top, modern_top$energy_pc, modern_top$`top_1_50%`)
# Archeological
#############################################################
# energy bounds
wd = gsub("Analysis", "Empirical Data/Energy Adaptation", dir)
setwd(wd)
bounds = fread("energy_boundaries.csv")
wd = gsub("Analysis", "Empirical Data/Kohler", dir)
setwd(wd)
kohler = fread("kohler.csv")
# merge with energy bounds
kohler_energy = merge(bounds, kohler, by = "Adaptation")
kohler_energy$beta = mapply(fit_beta, kohler_energy$`energy_50%`, kohler_energy$Gini)
# Pre-industrial
############################################################
wd = gsub("Analysis", "Empirical Data/Pre Industrial", dir)
setwd(wd)
pre_ind = fread("pre_industrial_ineq_energy.csv")
pre_ind$beta = mapply(fit_beta, pre_ind$energy.fit, pre_ind$Gini/100)
View(kohler)
merge = data.table(
source = c(modern$code, kohler_energy$Site, pre_ind$Country),
energy = c(modern$energy, kohler_energy$`energy_50%`, pre_ind$energy.fit),
beta = c(modern$beta, kohler_energy$beta, pre_ind$beta),
gini = c(modern$`gini_50%`, kohler_energy$Gini, pre_ind$Gini/100)
)
View(merge)
